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Abstract. We present a detailed investigation of neu- 
tron star atmospheres with low magnetic fields, B < 10® — 
10^° G, which do not affect opacities and equation of state 
of the atmospheric matter. We compute the atmospheric 
structure, emergent spectral fluxes and specific intensities 
for hydrogen, helium and iron atmospheres in a wide do- 
main of effective temperatures and gravitational acceler- 
ations expected for neutron stars. The iron atmospheres 
are computed with the opacities and equations of state 
from the OPAL opacity library. We show that the model 
atmosphere spectra are substantially different from the 
blackbody spectra. For light element atmospheres, the flux 
is greater than the blackbody flux, and the spectrum is 
harder, at high photon energies, whereas at low energies 
the spectral flux follows the Rayleigh- Jeans law with a 
(surface) temperature lower than the effective tempera- 
ture. The spectra of iron atmospheres display prominent 
spectral features in the soft X-ray range. The emergent 
specific intensity is anisotropic, with the anisotropy de- 
pending on energy. These properties of the atmospheric 
radiation should be taken into account for the proper in- 
terpretation of the thermal component of the neutron star 
radiation detectable at X-ray through UV energies. In par- 
ticular, fitting of the thermal component with the black- 
body model may result in substantially incorrect parame- 
ters of the radiating regions. We discuss an application of 
our atmosphere models to the nearby millisecond pulsar 
PSR J0437-4715. 

Key words: stars: neutron - dense matter - magnetic 
fields - radiative transfer - radiation mechanisms: thermal 



1. Introduction 

Thermal radiation from the surfaces of neutron stars 
(NSs) provides important information about their ther- 
mal histories and the properties of the high-density mat- 
ter in then interiors (Pethick 1992; Tsuruta 1995). The 
NS thermal-like radiation has been observed recently in 
the soft X-ray range with ROSAT (Ogelman 1995; and 
references therein) and EUVE (e. g., Foster et al. 1996; 
Halpern et al. 1996) and in the UV-optical range with 
HST (Pavlov et al. 1996). For the proper interpretation 
of these observations, one should take into account that 
NSs are not perfect blackbody emitters — properties of 
their radiation are determined, like in usual stars, by the 
radiative transfer in their atmospheres, and the emergent 
spectra can differ substantially from the Planck spectrum, 
depending on the chemical composition, magnetic field, 
energy flux and gravity in surface layers (Pavlov et al. 
1995, and references therein). 

There are at least two factors which make NS atmo- 
spheres different from those of normal stars: the immense 
gravity and huge magnetic flelds. The first factor is im- 
portant for all NSs; it leads to strong compression and 
non-ideality of the atmospheric matter, which affects con- 
siderably equation of state (EOS) and opacity. Even much 
greater and more complicated effects on the opacity and 
EOS can be expected from the very high magnetic field 
which, in particular, makes the NS atmospheres essentially 
anisotropic. 

First NS atmosphere models have been developed by 
Romani (1987; hereafter R87). He employed (nonmag- 
netic) opacities from the Los Alamos Opacity Library 
(LAOL) to calculate atmosphere models and emergent 
spectra for a few effective temperatures, 5.5 < log Tcff < 
6.5, and chemical compositions (pure helium, carbon, oxy- 
gen, iron, and a mixture of elements with cosmic abun- 
dance) and showed that the spectrum could indeed be 
quite different from the blackbody, especially for the he- 
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lium models. Further work has been concentrated on NS 
atmospheres with "typical" (strong) NS magnetic fields, 
B = 10^^ — 10^"^ G, which drastically change opacities and 
EOS (Miller 1992; Shibanov et al. 1992; Shibanov et al. 
1993; Ventura et al. 1993; Pavlov et al. 1994; Pavlov et 
al. 1995; Shibanov & Zavlin 1995; Shibanov et al. 1995a; 
Zavlin et al. 1995a). These papers (except that of Miller) 
were devoted to pure hydrogen NS atmospheres because 
no reliable opacities and EOS had been developed for non- 
ideal partly ionized gases of heavier elements. It has been 
demonstrated that the strong magnetic fields not only dra- 
matically affect the emergent spectra, but they also lead 
to remarkable anisotropy and strong polarization of the 
NS thermal radiation. 

Importance of the magnetic field depends on how large 
is the electron cyclotron energy, hfB = heB/nieC, in com- 
parison with other characteristic energies. For instance, 
if we consider a completely ionized (hot) atmosphere, 
where the main source of opacity is the free- free transitions 
and/or scattering on electrons, the opacity is not affected 
by the field for hvB <C hv, ksT (see, e. g., Kaminker et 
al. 1983). Since the energy transfer in the atmosphere is 
mainly driven by photons with hv ^ {\ — 10)fcBT, one 
can neglect the magnetic effects if i? ^ {mec/he)kBT ~ 
10^°{T/10^ K) G. If a substantial fraction of bound atoms 
or ions is present, the magnetic effects are determined by 
the ratio hvB/Z^Ry ~ (B/IO^ G)Z-'^, where Z is a char- 
acteristic ion charge, and Ry = 13.6 eV is the Rydberg 
energy. When this ratio is < 1, magnetic effects on the 
bound-bound transitions between lower quantum levels, 
as well as on the bound-free transitions well above pho- 
toionization thresholds, can be neglected. Although tran- 
sitions which involve highly excited levels may be affected 
by much lower fields, they usually are not important for 
conditions typical for NS atmospheres. 

Thus, depending on atmospheric temperature and 
chemical composition, one can neglect the magnetic effects 
on opacities and EOS if B < 10*^ - 10^° G. This means 
that there exists at least one class of NSs, millisecond pul- 
sars, for which the magnetic effects can be expected to 
be small. Four millisecond pulsars have been detected in 
soft X-rays (Halpern 1996, and references therein). The 
radiation from at least one of them - a nearby, bright 
PSR J0437 - 4715, with the field B - 3 x 10^ G - has 
been shown to have a thermal-like component apparently 
originating from the NS surface (Becker & Triimper 1993; 
Edelstein et al. 1995; Halpern et al. 1996). To interpret 
these data properly, the observed spectra and the light 
curves should be compared to nonmagnetic atmosphere 
models. One cannot exclude also that high magnetic fields 
of ordinary pulsars are concentrated in small "platelets" , 
whereas a substantial fraction of the NS surface has much 
lower fields (Page et al. 1995). Finally, since the NS mag- 
netic fields are subject to decay (although characteristic 
time of this decay remains highly uncertain — see, e. g., 
Bhattachrya 1995), one can expect that very old NSs, most 



of which do not show pulsar activity, have magnetic fields 
sufRciently low for nonmagnetic atmosphere models could 
be applied to analyze their thermal radiation. Detection of 
two old NS candidates has been reported recently (Stocke 
et al. 1995; Walter et al. 1996). Although identification 
of these objects as NSs requires further confirmation, the 
amount of dead pulsars observable through their thermal 
radiation should exceed considerably that of active pul- 
sars, so that one can expect more such objects to be dis- 
covered in near future. 

Thus, investigation of the nonmagnetic NS atmosphere 
models can be justified not only by their relative simplicity 
— the models can be applied to analysis of observations of 
real objects. Although a set of such models has been pre- 
sented in R87, there exist several reasons to re-examine 
the previous results and to continue work on the nonmag- 
netic models. 

(i) New, much improved opacity codes have been de- 
veloped by the OPAL group, and a comprehensive set of 
the opacity data has been computed (Iglesias et al. 1992; 
Rogers & Iglesias 1994). These new opacities substantially 
exceed those from LAOL for heavy elements just at the 
densities and temperatures characteristic for the NS at- 
mospheres. Particularly interesting for the application to 
NS atmospheres is a set of monochromatic opacities for 
pure iron composition, generated at request of, and in 
collaboration with, Roger Romani (Romani 1996, private 
communication). Rajagopal and Romani (1996; hereafter 
RR96; see also Romani et al. 1996) have used the results 
to develop modified NS atmosphere models. EJ 

(ii) No hydrogen atmosphere models and only three 
hchum models (for T^g = 10^■^ 10^ and lO'^-^) were 
presented in R87, and a few hydrogen and helium non- 
magnetic models were presented in our previous pa- 
pers (Shibanov et al. 1992; Zavlin et al. 1994; Pavlov et 
al. 1995), mainly for comparison with the magnetic mod- 
els, without a systematic analysis. However, the low-field 
NSs, including the millisecond pulsars, are assumed to 
have been passed, during their long life, through an ac- 
cretion stage when a substantial amount of light elements 
have been accumulated at the NS surface. (It may even 
be that the surface magnetic fields are low because of 
accretion- induced field decay — see, e. g., Bhattacharya 
1995). Due to gravitational stratification (Alcock & Illar- 
ionov 1980) and/or nuclear spallation reactions which de- 
stroy heavy elements (Bildstein et al. 1992, 1993), the up- 
per atmospheric layers can be expected to consist mainly 
of lightest elements unless they have been dragged down- 
ward by convection or burned in nuclear reactions. Only 
two hydrogen atmosphere models (for Tcff — 10^'^ and 
10^ K) were presented in RR96, which is not enough for 
a detailed analysis of how the atmospheric structure and 
emergent spectra depend on the model parameters. 



We became aware about RR96 after the present paper had 
been generally completed. 
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(iii) Only the (local) flux spectra of the emergent radi- 
ation have been studied in both R87 and RR96. However, 
the temperature distribution over the NS surface should 
be, as a rule, nonuniform due to either the presence of a 
hot polar cap in active pulsars (which is demonstrated 
by the light curve of PSR J0437 - 4715 — Becker & 
Triimper 1993) or anisotropy of the heat conduction in 
the NS crust (Shibanov & Yakovlev 1996). This means 
that the flux from a NS is to be determined by integra- 
tion of the angle-dependent specific intensities over the 
observable NS surface; this flux may differ substantially 
from the local flux due to the frequency-dependent limb- 
darkening effect which has not been investigated for the 
nonmagnetic models. In addition, gravitational bending 
of the photon trajectories is to be taken into account to 
obtain correct light curves of the NS thermal radiation 
(Shibanov et al. 1995b; Zavhn et al. 1995b, c). 

(iv) It would be desirable to verify the models pre- 
sented in R87 and RR96 using an alternative, more mod- 
ern numerical approach which has been applied by our 
group to modeling magnetic NS atmospheres. In partic- 
ular, that approach takes into account anisotropy of ra- 
diation in both solving the radiative transfer and calcu- 
lating the atmosphere temperature distribution, whereas 
the Lucy-Unsold method employed in R87 and RR96 for 
calculating the temperature distribution is correct only in 
the case of isotropic radiation field. In addition, Thomson 
scattering, treated in R87 with an oversimplified approach, 
and excluded from consideration in RR96, may be impor- 
tant for light element atmospheres at Tcs }t 10® K and is 
to be included consistently in the atmosphere modeling. 

In this paper we present a set of nonmagnetic NS at- 
mosphere models of different chemical compositions com- 
puted with an advanced atmosphere modeling code com- 
bined with the OPAL opacities. We study the atmospheric 
structure, and the spectra and angular distribution of 
emergent radiation in a wide range of effective temper- 
atures and surface gravities. We consider the atmospheres 
in radiative equilibrium, which means that the total en- 
ergy flux is transferred only by radiation, whereas other 
heat transport mechanisms (convection, electron heat con- 
duction) are of no importance. The models with convec- 
tion, which plays a role at lower Tcs, will be analyzed 
elsewhere (see Zavlin et al. 1996 for preliminary results). 
Our estimates show that, for the models considered, con- 
tribution of electron heat conduction to the heat flux never 
exceeds a few tenths of percent at unit optical depth. It 
may reach several percent only at very bottoms of heavy 
element atmospheres with lowest effective temperatures; 
the intensities of far Wien tails of the emergent radiation 
formed at these layers are too low to be of any practical 
importance. It should be mentioned that our estimates 
of the conductive opacity exceed those of RR96 (making 
electron heat conduction less important) because RR96 
used an approximate equation with the Coulomb loga- 
rithm omitted. On the other hand, the equations we use 



(Hochstim & Massel 1969) also become inaccurate at high 
densities and low temperatures, so that the problem needs 
further consideration. 

We describe our approach to atmosphere modeling in 
Section 2, and present the results in Section 3. Some ob- 
servational implications are discussed in Section 4. 

2. General description of the model 

2.1. Basic equations 

Since a characteristic thickness of the NS atmosphere, 
~ 0.01 — 1 cm, is much smaller than the NS radius, and 
characteristic densities are rather high, we consider the 
plane-parallel atmospheres in local thermodynamic equi- 
librium. Following a standard approach (e. g., Mihalas 
1978), we proceed from the radiative transfer equation in 
the form of a second-order boundary problem: 

fJ^^-^JI'-^Uiy ^ K{uu - S,,) , (1) 
ay Ki/ dy 

T-^-'^i^ = U^\y=0, Ui, = B„\y^oo , (2) 

Ki/ Ay 

where < /i < 1 is cosine of the angle between 
the normal to the atmosphere and the wave vector, 
Uu = u„{y,fi) = 5(/+ + /^), = I^{y,+n) and 

/j7 = Iu{y, "M) are the specific spectral intensities of the 
outward- and inward-directed streams of radiation, y is 
the column density (dy — —pdz), and B^, = B^{T(y)) 
is the Planck function. The monochromatic opacity k^, = 
k^(y) — a^{y) + cr^{y) is the sum of the true absorption 
opacity and scattering opacity a^. For the source func- 
tion Si, — Sy{y) we use 

Sv = {(^uJu + a^Bi,)k~^ , (3) 

where Jiy = Juiy) = /g^ Ui, d/i is the mean intensity. This 
form of the source function implies the isotropic scatter- 
ing approximation which proved to be very accurate even 
in the case when scattering dominates over true absorp- 
tion (e. g., Mihalas 1978). The boundary conditions (||) 
assume that there is no incident radiation at the outer 
boundary y = 0, and the photon energy density tends to 
its equilibrium value at large depths. 

In the present paper we consider only atmospheres 
in radiative equilibrium, where convection and electron 
heat conduction are negligible, and the energy is trans- 
ferred solely by radiation. In this case the radiative flux 
is independent of y and coincides with the total flux, 

/ F.d^^4vr/ dW df^f^^^asBT^s, (4) 
Jo Ja Jo K dy 

where F,j{y) is the spectral flux, usb is the Stefan- 
Boltzmann constant, and Tcff is the effective temperature. 
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It follows from Eqs. (|l]) and (Q) that the radiative equilib- 
rium condition can be also presented in the form 



Computing the radiation field starts fi'om calculating 
the Eddington factors, 



(5) /. = 



We consider the atmospheres in hydrostatic equilib- 
rium, which implies that the gas pressure is P = gy, where 
g is the gravitational acceleration (the radiative force is in- 
significant for the temperatures of interest, Tcs -C 10^ K). 
The above equations are to be supplemented with EOS, 
P = P(p,T), and true absorption and scattering opac- 
ities, a^{p,T) and (j^{p,T). The opacities can be either 
extracted from opacity libraries or calculated with the 
use of atomic data (radiative cross sections) . In the latter 
case, the set of equations of ionization equilibrium (EOIE) 
should be solved to find the fractions of ions in different 
stages of ionization as well as the electron number density 

Ne[p,T). 

The input parameters of a model are the effective tem- 
perature Toff, surface gravitational acceleration g (or the 
NS mass and radius) and chemical composition. Solving 
these equations yields the atmospheric structure {T{y), 
piy), iVe(y), etc.), and the radiative field in the atmo- 
sphere (F^{y), Jv(y), u^{y,p,), including the angular and 
spectral distribution of the emergent radiation, l^{0,p) — 
2u,{0,p)). 

It should be noted that the models give the radiation 
and the atmospheric structure for a local element of the 
NS surface. If the surface conditions are not uniform (e. g, 
the effective temperature may be higher at the magnetic 
poles), the observable photon flux should be obtained from 
integration of Ii,{Q,p) over the visible NS surface, with 
allowance for the gravitational redshift and bending of the 
photon trajectories (e. g., Shibanov et al. 1995b; Zavlin et 
al. 1995b, c). 

2. 2. Method of solution 

Solving the system of the model equations splits into two 
parts, each of which can be solved with the knowledge 
of the other: calculation of the radiation held given the 
structure of the atmosphere, and calculation the structure 
given the radiation field. We start from calculating the 
temperature distribution in the Eddington approximation, 
integrating the equation 



d T _ 3 
dy Teff 16 



kR 



Tcff 



(6) 



with the boundary condition r(0) ~ 0.841 Tgg. Here kfi — 
kji{T, p) is the Rosseland mean opacity which is calculated 
at each step of the integration, with p{y) being found from 
EOS and hydrostatic equilibrium. The initial atmospheric 
structure is used to obtain a first approximation for the 
radiation field. 



p Uu dp and h^, — 



pui, dp . (7) 



In deep layers, where the radiation is nearly isotropic and 
the diffusion approximation is valid, we have ~ ^ and 
hn 2± i. The dependences of fy and hy on y account for 
the depth dependence of anisotropy of radiation. Integrat- 
ing Eqs. (0) and ^ over p yields the ordinary differential 
equation for the mean intensity. 



1 T 1 '-^^^ {Ju Pu ) 

dy ky dy 

with the boundary conditions 



1 d 
ky dy 



(8) 



(9) 



Now fy{y) and hy{y) can be calculated by solving Eqs. (|]), 
(||) and (0) iteratively: we start from solving Eq. (|l|) with 
the source function = B^ to obtain an initial approxi- 
mation for fy and from Eq. (^; then we use these fi, 
and hv to solve (^), substitute the solution into the 
source function, and solve Eq. (|l]) to find new Ui, and, 
hence, fy and h^, and so on, until convergence of fy and 
hy. Note that at each step of these "inner" iterations 
we solve ordinary differential equations instead of origi- 
nal integro-differential (with respect to p and y) Eq. ([!]). 

In principle, one could use Jviy) obtained in the course 
of the inner iterations to calculate a new atmospheric 
structure. However, convergence of the "global" iterations 
can be significantly improved (Auer & Mihalas 1968) if in- 
stead we use the solution of the integro-differential (with 
respect to v and y) equation 



_d_J__d_ 

dy ky dy 

dPy Uyl (J, 



-a 



fvJv — Oil/ i^Ji/ Pu^ 
Py,) dv' 



(10) 



" dT ay,{dPy,/dT)dv' 



with the boundary conditions (||) and with the Eddington 
factors calculated as described above. Since the spectral 
ffux is 



Fy 



A-K d 
ku dy 



(11) 



the additional (as compared to Eq. (^) integral term in 
Eq. (10) automatically provides the total radiative flux F 
to be constant throughout the atmosphere at any given it- 
eration. Moreover, this integral term turns into zero iden- 
tically when the radiation field and atmospheric structure 
satisfy the radiative equilibrium condition (^). 

The mean intensity found from Eq. (10) is used to cal- 
culate a new temperature T + 5T for each atmospheric 
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layer, with the correction ST equal to the ratio of the two 
integrals in the additional term in Eq. (10). Once the new 
T{y) is obtained, we can complete calculation of the at- 
mospheric structure (and new opacities) making use of the 
hydrostatic equilibrium, EOS and EOIE, and then com- 
pute the corresponding radiation field as described above 
(starting with calculation of the Eddington factors). The 
global iterations continue until the maximum relative cor- 
rection \ST/T\ becomes smaller than a given value (typi- 
cally ^ 1%). Note that using the modified transfer equa- 
tion (10) substantially accelerates the global iterations (a 
reasonable accuracy is reached in 2-4 iterations). 

We solve the model equations numerically by discretiz- 
ing on the depth variable y, photon energy E = hv and 
angular variable ^. Typically, we use a grid with 100 — 150 
depth points logarithmically spaced between lO^'* — 10^^ 
and 1 — 10 g cm~^, 100 — 200 energy points logarithmi- 
cally spaced between 0.1 fcs T^s and 60 fcs Toff, and 50 an- 
gular points for /Lt between and 1. A fourth-order Hermi- 
tian discretization method (Auer 1976) is applied to solve 
Eqs. (H), (^ and (10). We also use the ehmination scheme 
suggested by Rybicki & Hummer (1991), which provides 
a higher accuracy in comparison with the standard Gaus- 
sian method. 

Since the monochromatic opacities of heavy elements 
(particularly iron) include many thousands of spectral 
lines, their frequency dependence should be properly 
binned to use them for construction of the atmosphere 
models. We average the original monochromatic opacities 
arithmetically in each bin of the logarithmic frequency 
grid (typically (1 — 2) x 10'^ frequency points per bin) and 
use these binned opacities to compute the atmospheric 
structure and a "smoothed" spectrum as described above. 
We checked by varying the number of bins that ^ 100—200 
bins provide sufficient accuracy (a relative error < 1 — 2%) 
of the atmospheric structure and the total flux. To obtain 
the spectrum with a fine spectral resolution, we calcu- 
late the intensities of the outward- and inward-directed 
streams, 

Ii^{y,fJ') = - [ S,yk^exp(-- f K dy2\ dyi (12) 



2.3. Opacities and EOS 

The true absorption opacity ai, is due to the free-free, 
bound-free and bound-bound transitions; CTj^ is dominated 
by the Thomson scattering on electrons. The main distinc- 
tion of the opacities, EOS and EOIE in a NS atmosphere 
from those in usual stars is due to high non-ideality (pres- 
sure effects). For modeling hydrogen/helium atmospheres, 
we calculate the opacities and EOIE of non-ideal gases 
using the occupation probability formalism (Hummer & 
Mihalas 1988) as described by Zavlin et al. (1994). The re- 
sults of these calculation agree fairly well with the OPAL 
data (Iglesias et al. 1992) obtained with a quite different 
approach (Rogers 1986). 

For modeling iron atmospheres, we use the opacities 
and EOS computed by the OPAL group, in collaboration 
with R. Romani. The opacities were computed for the 
domain logi? -5 (1) 1, logT = 4.5 (0.25) 7.5 (i? = p/T^, 
p in g cm~^, Tg = T/(10^ K)). They are provided on a 
grid of 10"* photon energies linearly spaced in the range 
2 X IQ-^ksT <E< 20kBT. To obtain k^ at given T and p, 
we use the cubic spline interpolation of log k^ with respect 
to logi? and logT. For E > 20kBT we use additional 
tables for the same R and T and 1 eV < i? < 10 keV. 



3. Results 

We have computed a set of the atmosphere models for 
NSs with various effective temperatures Toff, gravitational 
accelerations g and chemical compositions. In this paper 
we demonstrate representative models for pure hydrogen, 
helium and iron compositions. The hydrogen and helium 
models are presented for log Toff = 4.7 (0.3) 6.5; higher ef- 
fective temperatures are not expected to be observed in 
cooling NSs, and at lower Toff the helium atmospheres be- 
come convectively unstable. The iron atmosphere models, 
where convection can develop at higher Toff, are presented 
for log Toff = 5.3 (0.3) 6.5. We chose a "basic" gravita- 
tional accelerations 1714 = g/{10^'^ cm s~^) = 2.43 (it cor- 
responds to standard NS mass M — 1.4Mq and radius 
i? = 10 km); a few examples are given for 514 ~ 0.5 and 
4.5. 



and 



h = - / ^'^^'•^ exp (-- [ k^ dy2 ) dyi , (13) 

with the original monochromatic opacities (at ~ 10^ fre- 
quency points). The spectral flux is then computed as 



F4y) = 27r / [I+{y, m) ~ I'iy, ^i)]^,d^l . (14) 
"'0 

These "output" spectra can be further smoothed, if 
needed, with a standard spline procedure. 



3.1. Atmospheric structure 

The atmospheric structure has not been analizcd thor- 
oughly in the previous works; however, it is not only in- 
teresting in itself — it also allows one to better under- 
stand properties of the emergent radiation. Quite differ- 
ent structures of atmospheres composed of light and heavy 
elements are due to very different opacities. An example 
of the energy dependence of the opacities at the depths 
log y = —6, —4 and —2 is shown in Fig. 1 for the hydrogen, 
helium and iron atmosphere models with log Toff = 5.9, 
gii = 2.43. We see that the Fe opacities, dominated by 
numerous spectral lines, are, on average, greater than H 
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Fig. 1. Monochromatic opacities vs. photon energy at the layers y 
iron atmospheres with T^g = 10^'® K and g — 2.43 x 10^* cm s~^. The upper right panel shows the original Fe opacities (10^ 
energy points); the binned opacities (200 points) are shown in the bottom right panel. The horizontal dashed lines demonstrate 
the Rosseland mean opacities kn at the corresponding layers. 



and He opacities, and their mean slope at the most im- 
portant frequency range is almost flat. On the contrary, 
H and He opacities are dominated by the free-free and 
bound- free transitions {k^, oc oc E~^) at moderate 
and low frequencies (energies), whereas at high frequen- 
cies they are dominated by Thomson scattering (indepen- 
dent of ly) . These distinctions are responsible for different 
atmospheric structures and properties of the emergent ra- 
diation. 

Particularly important component of the atmospheric 
structure is the depth dependence of temperature, T{y). 
Examples of T{y) for several values of T^s and for three 
chemical compositions are demonstrated in Fig. 2. We see 
that, for any Tes, the temperature at a given depth y 
is systematically higher for iron atmospheres. In other 
words, a characteristic depth, j/off, which can be defined 
by the equation T(ycff) = T^s, is smaller for iron atmo- 
spheres (—4.6 < logyoff < —2.1 in our example) than for 
helium and hydrogen atmospheres (—3.7 < logj/off < 
and —2.9 < logj/cff < 0, respectively, for the same range 
of Toff). Higher temperatures for iron atmospheres at 
large depths, where T*{y) ~ T^sTr = T^g k^dy, are 
due to higher Rosseland mean opacity k^ (cf. Fig. 1). 
Substantially higher surface temperatures {Ts[= T(0)] = 
(0.78 - 0.81) Tefi for Fe vs. T, = (0.31 - 0.53) Teff for H, 
at log Teff = 5.3 — 6.5) can be explained by generally flat- 
ter frequency dependence of the Fe opacity, which results 
in temperature profiles and surface temperatures close to 
those for the grey atmosphere {Tg ~ 0.811 Toff). The dif- 
ferences between the hydrogen and helium temperature 
profiles depend on the degree of ionization in the atmo- 
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Fig. 2. Temperature profiles T{y) of the models with different 
effective temperatures and chemical compositions. The filled 
circles mark the points where T{y) = Tcs; the corresponding 
y = yeS can be considered as a characteristic depth of the 
atmosphere. 



sphere, which, in turn, is determined by Tofi and g. In par- 
ticular, at high Teff the light element plasma is completely 
ionized, and its opacity in deep, dense layers is mainly de- 
termined by the free-free transitions: k^ ~ Z^A~'^pT~'^/'^, 
where Z is the ion charge, and A the atomic weight. 
Using this expression for k^, Eq. (^, and EOS for the 
ideal gas, p oc A{Z + ly^gyT''^, we obtain T^^/MT oc 



V. E. Zavlin, G. G. Pavlov and Yu. A. Shibanov: Model neutron star atmospheres with low magnetic fields 



7 



-2 



CIO 

o 



1 1 1 1 1 


III 


1 1 1 1 j.-' 


- gi4=2.43 




/ / 
/ / - 






A"" 




/ / /// 

A 




/ // // 

/ // x// 

_ / // A 

- A'' A ' 

n 1 1 1 1 1 


' 

A 

III 


H - 

He - 

Fe — 

1 1 1 1 r 



-8 -6 -4 -2 
Log y [g cm"^] 

Fig. 3. Mass density profiles p(y) for two values of 
log Teff = 5.3 and 6.5 and three chemical compositions. The 
filled circles indicate the characteristic densities peft = pivea) 
(cf. Fig. 2). 



gZ^\A[Z + l)]^^?/dy, which yields the following equation 
for the asymptotic behavior of the temperature of a fully 
ionized atmosphere at large y: 



T{y) oc 



A 



\A 



. y4 



2/17 



V 



4/17 



(15) 



The value of the multiplier in front of y*^^^ for He ex- 
ceed that for H by 3% only; this means that the temper- 
ature profiles are essentially the same at sufficiently large 
depths, which we observe for the models with log Teff = 
6.5. At log Teff = 5.9, the hydrogen atmosphere is com- 
pletely ionized at the most important depths, whereas the 
helium atmosphere contains a few percent of once-ionized 
helium, Hell (see Fig. 4), which substantially increases the 
bound-free opacity (Fig. 1). This leads to a greater differ- 
ence between the He and H temperatures at y ~ ?/cff- At 
log Teff = 5.3 the difference is even more pronounced be- 
cause a larger fraction of HeH, > 10%, is available. Finally, 
at log Teff = 4.7, the H atmosphere contains ~ 10% of neu- 
tral atoms HI, and the He atmosphere is mainly composed 
of Hell, with an admixture of Hel. In both cases opacity 
is dominated by the bound-free transitions, and the tem- 
perature profiles are close to each other at j/ ~ j/off- 

The temperature dependence of the He atmosphere 
with log Teff = 4.7 shows an interesting anomaly at small 
y: it slightly increases outward at logy < — 7.5. This in- 
verse behavior is due to substantial growth with y of 
the Hel fraction at logy < — 7 (see Fig. 4), accompa- 
nied by increasing role of the Hel photoionization absorp- 
tion at > 20 — 25 eV. To provide the radiative equilib- 
rium (Eq. (H)) at small y, the local blackbody intensity 
B,y{T{y)) is to be greater than the mean intensity Ju{y) 
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Fig. 4. Depth dependences of ionization for models with dif- 
ferent values of log Tcff (numbers near the curves) . Two upper 
panels show profiles of non-ionized hydrogen and helium, xo{y), 
and once-ionized helium, a;i(y), for pure hydrogen and helium 
models. The bottom panel shows the depth dependences of the 
mean ion charge Z*{y) for pure iron models. The filled circles 
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at low energies (below the main photoionization jump), 
and smaller than Jv{y) at high E (above the jump). At 
logy 5; ~ 8 the main jump (provided by Hell) arises at 
~ 50 — 55 eV, whereas at —7 < logy < — 6 the jump is 
shifted to ~ 20 — 25 eV due to a greater fraction of Hel. 
Therefore, for energies between approximately 20 — 55 eV, 
Bi, > Ji/ at logy < — 8, while the opposite inequality holds 
at — 7 < logy < — 6. Since is practically independent 
of y in all these layers, the Planck function at logy < — 8 
should be higher than at — 7 < logy < — 6. This decrease 
of at these depths is possible if only the surface temper- 
ature T(0) is higher than the temperature in the deeper 
layers. 

Figure 3 shows the depth dependence of the mass 
density p(y) for the models with log T^s = 5.3 and 6.5. 
The density grows with y linearly (exponentially with the 
geometrical depth z) at the upper (exponential) atmo- 
sphere, where the temperature is close to Ts- In deep lay- 
ers, the density grows more gradually due to the steeper 
temperature growth. The density can be estimated as 
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p ^ tuhA^Z^ + l)^^gy{kBTY^, where Z, = Z^{y) = 
NeAmn / p is the mean ion charge defined as the ratio of 
the electron number density Ne to the summed number 
densities of atoms in all stages of ionization. Since the hy- 
drogen and helium atmospheres at the temperatures con- 
sidered are strongly ionized (Z* ~ Z), the density of a 
helium atmosphere at a given y exceeds that of a hydro- 
gen atmosphere with the same T^g and g. In particular, 
for deep layers of fully ionized atmospheres the curves p{y) 
has a universal slope, 

p OC + l)-15/17^19/17,^15/17y-^8/17yl3/17 (^g) 

SO that p^^'^\y) / p^^Hu) — 2-6. Iron atmospheres at these 
temperatures are only partially ionized, and the differ- 
ent disposition of the p^^'^\y) with respect to 
at log Tcff = 5.3 and 6.5 is explained by different de- 
grees of ionization (Z* ~ 5 — 7 and 21 — 22, respec- 
tively — see Fig. 4). Characteristic densities in different 
models, pcff = PiUes), are generally lower for iron at- 
mospheres, pcS ~ 0.003 — 0.02 g cm~'^, than for hydro- 
gen/hehum atmospheres, pcs ~ 0.005 — 1 g cm^'', for the 
range log Toff = 5.3-6.5. The ratio ZeS = ycs/peS gives 
a characteristic scale height in an atmosphere; it varies 
from ^ 0.01 cm for Fe at log Teff = 5.3 to 2 cm for H 
at log Tcft = 6.5. 

The atmospheric opacity and, hence, the temperature, 
density and properties of the emergent radiation depend 
essentially on the ionization state at different depths. Fig- 
ure 4 demonstrates the depth profiles of the fractions of 
ions in different stages of ionization, Xj = Nj/ {Nq -|- iVi + 
. . . -|- Nz) (j = 0, 1, . . . , Z is the ion charge). We show the 
nonionized fraction xq for the hydrogen models, nonion- 
ized and once-ionized fractions, xq and xi, for the helium 
models, and the mean ion charge Z» = J^j ^'^^ the 

iron models. The hydrogen atmospheres are strongly ion- 
ized at the considered temperatures (even at log Tes — 4.7 
the nonionized fraction does not exceed 17%); however, 
the contribution of the bound-free transitions into opacity 
is not negligible at xq as small as 10^'^ — 10^^ (cf. Fig. 1). 
At logy < — 2.5 the values of Xq are determined by the 
usual Saha equation: xq grows with y at very outer lay- 
ers, where the temperature is almost constant, due to 
increasing density; the growth decelerates in deeper lay- 
ers due to increasing temperature. At large depths, where 
P ^ 10^^ — 10^^ g cm^'^, the fraction xq decreases steeply 
due to pressure ionization. Similar behavior is seen for the 
helium fractions xq and xi , although the ionization degree 
is lower than for hydrogen due to higher ionization poten- 
tials. At moderate temperatures, the main contribution 
to the opacity comes from the bound-free transitions of 
Hell; the non- ionized Hel plays a role at lower Tcs- The 
iron atmospheres are only partly ionized even at high T^s 
and very deep layers (e. .g., Z* does not exceed 21 — 22 
through the whole atmosphere with log Tcff — 6.5) due to 
very high ionization potentials and binding energies. 



The atmospheric structure depends also on the gravi- 
tational acceleration g. Our computations show that vary- 
ing g in the range expected for the surfaces of NSs affects 
the temperature profile only slightly: T oc g^^^^ in very 
deep layers (see Eq. (p^), and the effect on the surface 
temperature is even smaller. The characteristic depth and 
density depend on g much stronger; in particular, p oc 5 in 
surface layers, and p cx g^^^^'^ in deep layers. The ioniza- 
tion degree is affected accordingly: generally, the ionized 
fractions grow with g, particularly at the surface layers. 

3. 2. Spectral flux 

In Fig. 5 we show the spectral fluxes of the emergent radi- 
ation for selected effective temperatures and three chemi- 
cal compositions. The blackbody fluxes TrB^{Tcg) are also 
shown for comparison. 

We see that the model spectra are substantially differ- 
ent from the blackbody spectra with the same T^s- The 
most striking distinction is much harder high-energy tails 
of the spectra emergent from H and He atmospheres (see 
also R87, Pavlov et al. 1995, and RR96). The reason for 
such behavior is that the H and He free-free and bound- 
free opacities rapidly decrease with increasing photon en- 
ergy (oc E^^ in a wide energy range — see Fig. 1), so 
that the radiation of higher energies is formed in deeper 
and hotter layers. The maxima of the spectra of strongly 
ionized atmospheres (at log Tcff > 5.0 and 5.5 for H and 
He, respectively) lie at E = (4.6 — 5.0)A:bT, i. e., they 
are shifted by a factor of 1.6 — 1.8 from the correspond- 
ing blackbody maxima towards higher energies. Thom- 
son scattering dominates the opacities at high energies 
(Fig. 1), which must be taken into account for calculating 
the high-energy tails (i. e., at E >1 keV for the models 
with log Teff = 5.9). 

At lower Tgg, the light element atmospheres are less 
ionized, and photoionization edges and spectral lines are 
seen in their spectra. In particular, the Lyman jump and 
Lyman-a line of Hell become perceptible in the helium 
atmosphere spectra at log Toft J; 6, when the fraction of 
the once-ionized helium is as small as a few percent (see 
Fig. 4). The same fraction of neutral hydrogen is needed 
to form the hydrogen Lyman features which become well 
pronounced at log Teff < 5. 

The low-energy tails [E -^ksT) oi the model spectra 
follow the Rayleigh- Jeans law (oc E"^) corresponding to the 
surface temperature Ts{< Toff), i. e., F^/irBi, = T^j Toff < 
1. This occurs because the opacities grow with decreasing 
E so that the low-energy radiation emerges from the very 
surface layers. The low-energy flux is most strongly sup- 
pressed (by a factor of up to ~ 3.2) in hot H and He 
atmospheres because of smallness of the ratio Tsj Toff (see 
Fig. 2). 

The spectra emitted from the iron atmospheres are 
much more complicated, mainly due to numerous spectral 
lines and photoionization edges produced by ions in vari- 
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ous stages of ionization. Generally, they are closer to the 
blackbody spectra, and their high-energy tails are, on av- 
erage, softer than those emitted from light element atmo- 
spheres; their low-energy tails lie only slightly below the 
blackbody (Rayleigh-Jeans) spectra. The reason is that 
the energy dependence of the iron opacities, in spite of 
its complexity, is, on average, flatter than that for H and 
He. This makes the emergent spectrum to be closer to the 
grey atmosphere spectrum which, in turn, only slightly 
differs from the blackbody spectrum. However, the local 
deviations from the blackbody spectra are quite substan- 
tial, especially in the regions of the photoionization edges, 
and strongly dependent on the effective temperature. 

Figure 6 demonstrates how the gravitational accelera- 
tion affects the spectral flux emitted by the helium atmo- 
sphere with log Tcff = 5.0. The effect is seen most vividly 
in the Lyman continuum of Hell {E > 54 eV): the Ly- 
man discontinuity becomes progressively smoothed, and 
the trough a± E 0.1 keV becomes shallower, with in- 
creasing g. This occurs because the density grows with 
g and the pressure effects smooth out the photoioniza- 
tion jump in the opacity spectra (Zavlin et al. 1994). The 
gravity effect (pressure broadening) is also seen in the pro- 
flies of the absorption lines, particularly of the 20 eV line 
of Hcl. At higher temperatures, when the light element 
atmospheres are strongly ionized and the spectra are con- 
tinuous, the gravity effect becomes insigniflcant. In the 
case of iron composition the effect of surface gravity can 
be seen at higher temperatures. 
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3.3. Specific intensities: spectra and angular distributions 

Figure 7 demonstrates spectra of the specific intensities of 
the emergent radiation at different angles 9 = arccos fi be- 
tween the direction of emission and the surface normal for 
the models with hydrogen and iron compositions and the 
effective temperature of 10^'^ K. We see that the inten- 
sity is rather anisotropic — it decreases with increasing 
9 (limb-darkening effect). This is explained by the fact 
that radiation emitted at larger 9 escapes from shallower 
atmospheric layers where the temperature is lower. The 
anisotropy depends on photon energy (the intensity spec- 
trum varies with 9) and on chemical composition of the 
atmosphere. Note that the maximum of the intensity spec- 
trum shifts towards lower energies with increasing 9. 

The energy-dependent degree of anisotropy can be de- 
scribed by the ratio 



/.(m^ l)-/.(/i = 0) 
= 1) + = 0) 



(17) 



plotted in Fig. 8; it varies from to 1 for atmospheres 
with the temperature growing inward monotonously. For 
the light element atmospheres, the anisotropy is small 
at low energies because the radiation emerges from sur- 
face layers at any 6. It grows with E because radia- 
tion escapes from hotter layers at = 0. At very large 
E, when opacity is dominated by Thomson scattering, 
the angular dependence of the intensity is determined by 
the Ambartsumian-Chandrasekhar function H(/i) (e. g., 
Chandrasekhar 1960), and the degree of anisotropy is 
Qi, = 0.488 (qualitatively, one can use the diffusion ap- 
proximation, which yields I^{fJ.) oc l + 1.5fi and ~ 0.43). 
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Fig. 8. Degree of anisotropy a^, (Eq. (|l7|)) vs. photon energy 
for radiation emergent from hydrogen and iron atmospheres 
with the same parameters as in Fig. 7. 



In the case of iron atmospheres, the anisotropy is quite 
different. In particular, it is very small in the regions of 
strong spectral lines (e. g., around 0.2 and 1 keV) which 
originate from the surface layers. (Similar reduction of 
anisotropy near the spectral lines takes place for the light 
element atmospheres at lower Tcff.) The radiation from 
the iron atmospheres is very anisotropic {a^, — *■ 1) at 
high energies because the opacity keeps decreasing with 
increasing E (Thomson scattering is negligible — see 
Fig. 1). 
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Figure 9 shows representative angular distributions of 
the specific intensity emergent from hydrogen, helium and 
iron atmospheres at different photon energies. Although 
these distributions do not show the prominent features 
characteristic of radiation emitted from strongly magne- 
tized NS atmospheres (see Fig. 4 in Shibanov & Zavlin 
1995), this picture clearly illustrates strong dependences 
of the distributions on both the energy and chemical com- 
position. 

4. Discussion 

The presented examples of the atmosphere models show 
that the atmospheric structure and the properties of 
emergent radiation depend substantially on the chemi- 
cal composition of radiating layers and on the value of 
the effective temperature Teff- For the light element at- 
mospheres with a relatively high effective temperatures 
(log Teff > 5.3 for hydrogen, log Teff > 6 for helium), the 



atmospheric matter is strongly ionized, the (continuous) 

spcctriim is much harder at X-ray energies than the black- 
body spectrum, and the flux is lower than the blackbody 
flux at low (UV/optical) energies. For lower Teff, spec- 
tral features arise in the emergent spectra, particularly 
in the UV/EUV domain. These distinctions mean that 
fitting thermal spectra observed from NSs with a black- 
body model may be quite misleading. In particular, the 
blackbody fit of the high-energy (X-ray) tail of an ob- 
served spectrum would yield systematically overestimated 
effective temperature and underestimated A/d^ ratio {A is 
the area of the radiating surface, and d is the distance). 
The error introduced depends on the energy range and 
response function of the specific detector; c. g., for the 
ROSAT PSPC the temperature may be overestimated by 
a factor of up to ~ 3, and A/di^ underestimated by more 
than an order of magnitude, at detectable NS tempera- 
tures. On the contrary, at low (UV-optical) energies T^g 
is underestimated or A/di^ overestimated, typically by a 
factor of 2 — 3. To obtain reliable estimates of these quan- 
tities, one should fit the results with various model atmo- 
sphere spectra. It may happen that several models (e. g., 
with different chemical compositions) would equally well 
fit the data in a narrow spectral range. To avoid the ambi- 
guity, multiwavclcngth (c. g., optical through X-ray) ob- 
servations should be carried out, whenever possible, and 
fitted with the atmosphere models. 

Spectra of heavy element atmospheres are much more 
complicated — mainly because of numerous spectral lines, 
as we see from the examples of iron atmospheres. Note 
that although the shapes of our model spectra are qualita- 
tively similar to those obtained in RR96, the quantitative 
differences are significant, especially in domains around 
the strong photoionization features. The difference may 
arise from both the different procedures for opacity bin- 
ning and different numerical approaches used for the mod- 
eling (e. g., due to the inaccuracy of the Lucy-Unsold tem- 
perature correction procedure used in RR96). Although 
the overall shape of the iron atmosphere spectra is not 
so drastically different from the blackbody, using black- 
body models for fitting in narrow energy intervals may be 
even more misleading than for light element atmospheres 
because even the smoothed shape of the atmosphere spec- 
trum depends very sharply on Toff. For instance, we see 
from Fig. 5 that increasing Toff by a factor of 2 completely 
changes the shape of the iron atmosphere spectrum, and 
fitting the spectra with the blackbody in a range of, e. g., 
0.2 — 0.6 keV would lead to quite different (and quite ir- 
relevant) results. 

We do not know a priori what is the chemical compo- 
sition of the NS surfaces. We can only expect that if there 
are light elements, they should concentrate at the very 
surface, but we cannot exclude that the light elements are 
virtually absent, and the surface layers are composed of 
heavier elements. To investigate the chemical composition, 
observations with higher spectral resolution in the soft X- 
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ray range would be very useful. Our examples of the iron 
models show that the main spectral features (groups of 
lines and/or photoionization edges) can be resolved with 
a quite moderate resolution, ~ 3 — 30, which can be eas- 
ily achieved with modern CCD detectors. Thus, soft X-ray 
observations of NSs with future missions {AXAF, XMM, 
ASTRO-E) and their interpretation with the NS atmo- 
sphere models will be able to determine directly which 
elements reside at the NS surface. 

Our results also show a very important property of the 
atmospheric radiation — it is anisotropic, with anisotropy 
depending on frequency. This dependence is different for 
different chemical compositions and effective tempera- 
tures. This means that if the temperature is not uniform 
over the NS surface (which is very plausible for real NSs), 
one not only should take the anisotropy into account to 
explain the light curves, but also one cannot fit properly 
the spectra without allowance for anisotropy. To calculate 
the flux as measured by a distant observer, one should 
integrate the specific fiuxes, I^cosd, over the visible NS 
surface with allowance for the gravitational redshift and 
bending the photon trajectories. For instance, if the radi- 
ation originates from a small hot spot (polar cap) with a 
uniform Tcff and chemical composition, the flux as mea- 
sured by a distant observer is given by the equation (Zavlin 
et al. 1995c) 

A 1 

d^(TTIF^''°^^'^' ^^^^ 

where Aa is the apparent area of the hot spot, z = [I — 
{2GM/c^R)]-^/'^ - 1 is the gravitational redshift, = 
(l+z)i^, /ic is the cosine of the angle 9c between the radius- 
vector of the spot center and wave vector of the photon 
which travels to the observer along the bended trajectory. 
For reahstic NS radii, R > 5.6(M/Mq) km {z < 0.45), the 
apparent area can be estimated as Aa — irR^j^f^c, where 
7 is the angular size of the spot. The important feature 
of Eq. (|^ is that the observed flux is proportional not to 
the local flux, but to the speciflc intensity whose spectrum 
(see Figs. 5 and 7) is substantially different from that of 
the redshifted local flux used in RR96 (their Eq. (7)). 

The low-field atmosphere models can be directly ap- 
plied to the nearby {d ~ 140 pc) millisecond pulsar PSR 
J0437 - 4715 (P = 5.75 ms) detected with the ROSAT 
PSPC (Becker & Triimper 1993) and EUVE DS (Edel- 
stein et al. 1995; Halpern et al. 1996). Fitting the ROSAT 
spectrum of this low-field (i? ~ 3 x 10^ G) pulsar with 
the single-component blackbody model proved to be un- 
satisfactory — the best-fit blackbody fiux (correspond- 
ing to T = 1.5 X 10^ K) is lower than the observed flux 
in high-energy PSPC channels (above ~ 1 keV). How- 
ever, the ROSAT spectrum can be satisfactorily fitted 
with double-component models. For the "power-law + 
blackbody" fit, Becker & Triimper (1993) obtained T = 
1.7 X 10^ K and the radiating area A — 0.05 km^, while 
Halpern et al. (1996) obtained T = (1.0 - 3.3) x 10^ K 



and A = 0.008 — 1.1 km^. In this model, the power law 
spectrum is of a nonthermal origin, whereas the blackbody 
spectrum is produced by the hot polar cap of the pulsar. It 
remains unclear with this interpretation why there is no a 
phase shift between the pulsations at low energies (power 
law component) and high energies (thermal component) 
— being generated at different distances from the NS sur- 
face, these components are expected to be phase-shifted 
and to have different shapes. The values of the emitting 
area are also uncomfortably low in comparison with that 
expected from standard models for the radio pulsar po- 
lar cap, A = 27r2i?3/(cP) - 10 km^ (Arons 1981). Edel- 
stein et al. (1995) ruled out that this two-component fit is 
compatible with the EUVE data, whereas Halpern et al. 
(1996) came to the opposite conclusion and explained the 
EUVE flux with the power law component. 

It is natural to expect that fitting of the same data 
with the model atmosphere spectra and light curves would 
yield quite different results. Indeed, our preliminary anal- 
ysis of both the ROSAT and EUVE data indicates that a 
single-component model, in which the soft X-ray radiation 
originates from two polar caps of areas ^4 ~ 2 — 3 km^ cov- 
ered with hydrogen or helium with Tcs ~ (0.9 — 1.0) x 10^ 
K, fits satisfactorily both the spectra and the light curves 
(Becker et al. 1996). These results differ substantially from 
those of RR96 who fitted the J0437-4715 spectrum with 
the local flux spectra, without allowance for the limb dark- 
ening and gravitational lensing. It worths noting that the 
atmosphere models can be used to predict the flux from 
the NS in the far-UV range (A < 2000 - 3000 A), where it 
should overtake the flux of the cold (T ~ 4000 K) white 
dwarf companion. Measuring this flux (e. g., with HST) 
would be a critical test to distinguish between different 
models of the radiating region. Another important test of 
the models would be to observe the pulsar radiation in 
the soft X-ray range with a spectral resolution sufficient 
to resolve spectral features; observations at moderate X- 
ray energies (~ 1 — 10 keV) would be useful to constrain 
the power law component. Such observations would be 
quite feasible with future X-ray missions {SAX, AXAF, 
XMM, ASTRO-E). 
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